library(tidyverse)
library(ggdist)
library(ggside)
library(easystats)
library(patchwork)
library(brms)

logmod <- function(x) sign(x) * log(1 + abs(x))
sqrtmod <- function(x) sign(x) * sqrt(abs(x))
cbrtmod <- function(x) sign(x) * (abs(x)**(1 / 3))

df <- read.csv("../data/preprocessed_illusion1.csv") |>
  mutate(
    Block = as.factor(Block),
    Illusion_Side = as.factor(Illusion_Side)
  ) 
for (ill in c(
  "Ebbinghaus",
  "VerticalHorizontal",
  "MullerLyer"
)) {
  print(ill)
  model <- brms::brm(
    Error ~ t2(Illusion_Difference, Illusion_Strength, bs = "tp") + (1 | Participant),
    data = filter(df, Illusion_Type == ill),
    family = "bernoulli",
    algorithm = "sampling",
    refresh = 400
  )

  name <- paste0("gam_", tolower(ill), "_err")
  assign(name, model) # Rename with string
  save(list = name, file = paste0("models/", name, ".Rdata"))

  model <- brms::brm(
    brms::bf(
      RT ~ t2(Illusion_Difference, Illusion_Strength, bs = "tp") + poly(ISI, 2) + (1 | Participant),
      sigma ~ t2(Illusion_Difference, Illusion_Strength, bs = "tp") + poly(ISI, 2) + (1 | Participant),
      beta ~ t2(Illusion_Difference, Illusion_Strength, bs = "tp") + poly(ISI, 2) + (1 | Participant)
    ),
    data = filter(df, Illusion_Type == ill, Error == 0),
    family = "exgaussian",
    init = 0,
    algorithm = "sampling",
    refresh = 400
  )

  name <- paste0("gam_", tolower(ill), "_rt")
  assign(name, model) # Rename with string
  save(list = name, file = paste0("models/", name, ".Rdata"))
}

Ebbinghaus

Model Selection

test_models <- function(data) {
  # TODO: add random effect
  models_err <- list()
  models_rt <- list()
  for (diff in c(
    "Illusion_Difference",
    "logmod(Illusion_Difference)",
    "sqrtmod(Illusion_Difference)",
    "cbrtmod(Illusion_Difference)"
  )) {
    for (ill in c(
      "abs(Illusion_Strength)",
      "logmod(abs(Illusion_Strength))",
      "sqrtmod(abs(Illusion_Strength))",
      "cbrtmod(abs(Illusion_Strength))"
    )) {
      print(paste(diff, "x", ill))
      err <- glmmTMB::glmmTMB(
        as.formula(
          paste0(
            "Error ~ Illusion_Effect / (",
            diff,
            " * ",
            ill,
            ") + (1|Participant)"
          )
        ),
        data = data,
        family = "binomial"
      )
      models_err[[paste(diff, "*", ill)]] <- err

      rt <- glmmTMB::glmmTMB(
        as.formula(
          paste0(
            "RT ~ Illusion_Effect / (",
            diff,
            " * ",
            ill,
            ") + poly(ISI, 2) + (1|Participant)"
          )
        ),
        data = data
      )
      models_rt[[paste(diff, "*", ill)]] <- rt
    }
  }
  mutate(performance::test_performance(models_err), Outcome = "Error") |>
    rbind(mutate(performance::test_performance(models_rt), Outcome = "RT")) |>
    select(-Model, -log_BF) |>
    datawizard::convert_na_to(select = "BF", replacement = 1) |>
    arrange(Outcome, desc(BF)) |>
    display(footer = "Each model is compared to 'Illusion_Difference'")
}

test_models(filter(df, Illusion_Type == "Ebbinghaus"))

Error Rate

data <- filter(df, Illusion_Type == "Ebbinghaus")

Descriptive

plot_desc_errors_obs <- function(data, n_groups=5) {
  data |>
    mutate(
      Illusion_Difference = datawizard::categorize(Illusion_Difference,
        n_groups = n_groups,
        split = "equal_length",
        label = "median"
      ),
      bins = datawizard::categorize(Illusion_Strength, split = "equal_length", n_groups = 20)
    ) |>
    group_by(Illusion_Difference, bins) |>
    summarize(
      Illusion_Strength = mean(Illusion_Strength),
      Error = sum(Error) / n()
    ) |>
    ggplot(aes(x = Illusion_Strength)) +
    geom_bar(aes(y = Error, fill = Illusion_Difference),
      alpha = 1 / 3, position = "dodge", stat = "identity", width = diff(range(data$Illusion_Strength)) / 20
    )
}


plot_desc_errors <- function(data) {

  dat <- data |>
    mutate(
      Illusion_Difference =
        datawizard::categorize(Illusion_Difference,
          n_groups = 5,
          split = "equal_length",
          label = "median"
        )
    )

  col <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(dat$Illusion_Difference)))

  plot_desc_errors_obs(data) +
    geom_smooth(data=dat,
                aes(y = Error, color = Illusion_Difference, fill = Illusion_Difference),
      method = "gam",
      formula = y ~ s(x, bs = "cr", k=8),
      method.args = list(family = "binomial"),
      alpha = 1 / 3
    ) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
    scale_color_manual(values = col) +
    scale_fill_manual(values = col) +
    coord_cartesian(ylim = c(0, 1), xlim = range(dat$Illusion_Strength)) +
    labs(x = "Illusion Strength", y = "Probability of Error") +
    theme_modern() +
    facet_grid(~Illusion_Side, labeller = "label_both")
}

plot_desc_errors(data)

Model Specification

formula <- brms::bf(
  Error ~ Illusion_Effect / (logmod(Illusion_Difference) * abs(Illusion_Strength)) +
    (1 + Illusion_Effect / (logmod(Illusion_Difference) * abs(Illusion_Strength)) | Participant),
  family = "bernoulli",
  decomp = "QR"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

illusion1_ebbinghaus_err <- brms::brm(formula,
  data = data,
  refresh = 100,
  normalize = FALSE
)

save(illusion1_ebbinghaus_err, file = "models/illusion1_ebbinghaus_err.Rdata")

Model Inspection

load("models/illusion1_ebbinghaus_err.Rdata")
load("models/gam_ebbinghaus_err.Rdata")
plot_model_errors <- function(data, model, gam) {
  pred <- estimate_relation(model, length = c(3, 100))

  dat <- data |>
    mutate(
      Illusion_Difference = datawizard::categorize(
        Illusion_Difference,
        n_groups = 3,
        split = "equal_length",
        label = "median"
      ),
      Illusion_Strength = as.character(Illusion_Strength)
    ) |>
    group_by(Illusion_Strength, Illusion_Difference) |>
    summarize(
      Error = sum(Error),
      n = n()
    ) |>
    group_by(Illusion_Strength) |>
    mutate(
      n = sum(n),
      Error = Error / n,
      Illusion_Strength = as.numeric(Illusion_Strength)
    ) |>
    ungroup()

  plot_desc_errors_obs(data, n_groups = 3) +
    scale_fill_manual(values = c("#F44336", "#FFC107", "#4CAF50"), guide = "none") +
    ggnewscale::new_scale_fill() +
    geom_ribbon(data=pred, aes(
      ymin = CI_low, ymax = CI_high,
      fill = Illusion_Difference, group = Illusion_Difference
    ),
    alpha = 0.2
    ) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    geom_hline(yintercept = c(0.5), linetype = "dotted", alpha = 0.5) +
    geom_line(
      data = estimate_relation(gam, length = c(3, 100)),
      aes(y=Predicted, color = Illusion_Difference, group = Illusion_Difference),
      linetype = "dotted"
    ) +
    geom_line(data=pred, aes(y=Predicted, color = Illusion_Difference, group = Illusion_Difference)) +
    scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_color_gradientn(colours = c("#4CAF50", "#FFC107", "#F44336"), trans = "reverse") +
    scale_fill_gradientn(colours = c("#4CAF50", "#FFC107", "#F44336"), trans = "reverse") +
    theme_modern(axis.title.space = 5) +
    guides(color = "none") +
    labs(
      color = "Difficulty",
      fill = "Difficulty",
      y = "Probability of Error",
      x = "Illusion Strength"
    ) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
}

p_ebbinghaus_err <- plot_model_errors(data, illusion1_ebbinghaus_err, gam_ebbinghaus_err)
p_ebbinghaus_err

Model Performance

performance::performance(illusion1_ebbinghaus_err, metrics = c("R2", "ICC")) |>
  display()
R2 R2 (marg.) ICC
0.28 0.16 0.06

Response Time

data <- filter(df, Illusion_Type == "Ebbinghaus", Error == 0)

Descriptive

plot_desc_rt <- function(data) {
  dat <- data |>
    mutate(
      Illusion_Difference =
        datawizard::categorize(Illusion_Difference,
          n_groups = 5,
          split = "equal_length",
          label = "median"
        )
    )
  
  col <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(dat$Illusion_Difference)))

  dat |>
    ggplot(aes(x = Illusion_Strength, y = RT)) +
    geom_jitter2(aes(color = Illusion_Difference),
      width = diff(range(dat$Illusion_Strength)) / 50,
      alpha = 1 / 3
    ) +
    geom_smooth(aes(y = RT, color = Illusion_Difference, fill = Illusion_Difference),
      method = "gam",
      formula = y ~ s(x, bs = "cr", k=8),
      alpha = 1 / 10
    ) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    scale_color_manual(values = col) +
    scale_fill_manual(values = col) +
    coord_cartesian(ylim = c(0, 2)) +
    labs(x = "Illusion Strength", y = "Response Time") +
    theme_modern() +
    ggside::geom_ysidedensity(aes(fill = Illusion_Difference), color = NA, alpha = 0.3) +
    ggside::theme_ggside_void() +
    ggside::scale_ysidex_continuous(expand = c(0, 0)) +
    ggside::ggside() +
    facet_grid(~Illusion_Side, labeller = "label_both")
}

plot_desc_rt(data)

Model Specification

formula <- brms::bf(
  RT ~ Illusion_Effect / (Illusion_Difference * abs(Illusion_Strength)) + poly(ISI, 2) +
    (1 + Illusion_Effect / (Illusion_Difference * abs(Illusion_Strength)) | Participant),
  sigma ~ Illusion_Effect / (Illusion_Difference * abs(Illusion_Strength)) + poly(ISI, 2) +
    (1 | Participant),
  beta ~ Illusion_Effect / (Illusion_Difference * abs(Illusion_Strength)) + poly(ISI, 2) +
    (1 | Participant),
  family = "exgaussian",
  decomp = "QR"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

illusion1_ebbinghaus_rt <- brms::brm(formula,
  data = data,
  refresh = 100,
  init = 0,
  normalize = FALSE
)

save(illusion1_ebbinghaus_rt, file = "models/illusion1_ebbinghaus_rt.Rdata")

Model Inspection

load("models/illusion1_ebbinghaus_rt.Rdata")
load("models/gam_ebbinghaus_rt.Rdata")
plot_model_rt <- function(data, model, gam) {
  pred <- estimate_relation(model, at = c("Illusion_Effect", "Illusion_Difference", "Illusion_Strength"), length = c(3, 100))

  dat <- data |>
    mutate(Illusion_Difference = datawizard::categorize(
      Illusion_Difference,
      n_groups = 3,
      split = "equal_length",
      label = "median"
    ))

  col <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(dat$Illusion_Difference)))

  pred |>
    ggplot(aes(x = Illusion_Strength, y = Predicted)) +
    ggdist::stat_gradientinterval(
      data = dat,
      aes(y = RT, fill = Illusion_Difference),
      geom = "slab",
      scale = 1,
      width = diff(range(dat$Illusion_Strength)) / 20,
      fill_type = "gradient",
      position = "dodge"
    ) +
    scale_fill_manual(values = col, guide = "none") +
    ggnewscale::new_scale_fill() +
    geom_ribbon(aes(
      ymin = CI_low, ymax = CI_high,
      fill = Illusion_Difference, group = Illusion_Difference
    ),
    alpha = 0.2
    ) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    geom_line(aes(color = Illusion_Difference, group = Illusion_Difference)) +
    geom_line(
      data = estimate_relation(gam, at=c("Illusion_Difference", "Illusion_Strength"), length = c(3, 100)),
      aes(color = Illusion_Difference, group = Illusion_Difference),
      linetype = "dotted"
    ) +
    scale_y_continuous(expand = c(0, 0)) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_color_gradientn(colours = rev(col), trans = "reverse") +
    scale_fill_gradientn(colours = rev(col), trans = "reverse") +
    coord_cartesian(ylim = c(0, 2.5)) +
    labs(
      x = "Illusion Strength",
      y = "Response Time (s)",
      color = "Difficulty",
      fill = "Difficulty"
    ) +
    theme_modern(axis.title.space = 5)
  # ggnewscale::new_scale_fill() +
  # scale_fill_manual(values = c("#F44336", "#FFC107", "#4CAF50"), guide = "none") +
  # ggside::geom_ysidedensity(data=dat,
  #                           aes(fill = Illusion_Difference, y = RT), color = NA, alpha = 0.3) +
  # ggside::theme_ggside_void() +
  # ggside::scale_ysidex_continuous(expand = c(0, 0)) +
  # ggside::ggside()
}

p_ebbinghaus_rt <- plot_model_rt(data, illusion1_ebbinghaus_rt, gam_ebbinghaus_rt)
p_ebbinghaus_rt

Model Performance

performance::performance(illusion1_ebbinghaus_rt, metrics = c("R2", "ICC")) |>
  display()
R2 R2 (marg.)
0.41 0.03
plot_ppcheck <- function(model, gam) {

  pred <- modelbased::estimate_prediction(model, keep_iterations = 50) |>
    bayestestR::reshape_iterations() |>
    mutate(iter_group = as.factor(iter_group)) |>
    estimate_density(select = "iter_value", at = "iter_group")

  predgam <- modelbased::estimate_prediction(gam, keep_iterations = 50) |>
    bayestestR::reshape_iterations() |>
    mutate(iter_group = as.factor(iter_group)) |>
    estimate_density(select = "iter_value", at = "iter_group")

  estimate_density(insight::get_data(model)$RT) |>
    ggplot(aes(x = x, y = y)) +
    geom_area(fill = "#9E9E9E") +
    geom_line(
      data = predgam,
      aes(group = iter_group), color = "#2196F3", size = 0.1, alpha = 0.5
    ) +
    geom_line(
      data = pred,
      aes(group = iter_group), color = "#FF5722", size = 0.1, alpha = 0.5
    ) +
    scale_y_continuous(expand = c(0, 0)) +
    scale_x_continuous(expand = c(0, 0)) +
    coord_cartesian(xlim = c(0, 2)) +
    theme_modern() +
    labs(x = "Reaction Time (ms)", y = "", title = "Posterior Predictive Check") +
    theme(
      plot.title = element_text(face = "bold", hjust = 0.5),
      axis.text.y = element_blank()
    )
}

plot_ppcheck(illusion1_ebbinghaus_rt, gam_ebbinghaus_rt)

Visualization

plot_all <- function(data, p_err, p_rt, question="") {
  illname <- unique(data$Illusion_Type)
  # Get stimuli
  dat <- data |>
    filter(Error == 0) |>
    filter(
      Illusion_Type == illname,
      Answer %in% c("left", "up")
    ) |>
    select(Stimulus, Illusion_Strength, Illusion_Difference)

  dat <- rbind(
    filter(dat, Illusion_Strength == min(Illusion_Strength)) |>
      filter(Illusion_Difference %in% c(min(Illusion_Difference), max(Illusion_Difference))),
    filter(dat, Illusion_Strength == max(Illusion_Strength)) |>
      filter(Illusion_Difference %in% c(min(Illusion_Difference), max(Illusion_Difference))),
    filter(dat, Illusion_Difference == min(Illusion_Difference)) |>
      filter(Illusion_Strength %in% c(min(Illusion_Strength), max(Illusion_Strength))),
    filter(dat, Illusion_Difference == max(Illusion_Difference)) |>
      filter(Illusion_Strength %in% c(min(Illusion_Strength), max(Illusion_Strength)))
  )


  img_leftdown <- filter(dat, Illusion_Difference == max(Illusion_Difference)) |>
    filter(Illusion_Strength == min(Illusion_Strength)) |>
    pull(Stimulus) |>
    unique()
  img_rightdown <- filter(dat, Illusion_Difference == max(Illusion_Difference)) |>
    filter(Illusion_Strength == max(Illusion_Strength)) |>
    pull(Stimulus) |>
    unique()
  img_leftup <- filter(dat, Illusion_Strength == min(Illusion_Strength)) |>
    filter(Illusion_Difference == min(Illusion_Difference)) |>
    pull(Stimulus) |>
    unique()
  img_rightup <- filter(dat, Illusion_Strength == max(Illusion_Strength)) |>
    filter(Illusion_Difference == min(Illusion_Difference)) |>
    pull(Stimulus) |>
    unique()


  img_leftdown <- paste0("../session1/stimuli/", img_leftdown, ".png") |>
    png::readPNG() |>
    grid::rasterGrob(interpolate = TRUE) |>
    patchwork::wrap_elements()
  img_rightdown <- paste0("../session1/stimuli/", img_rightdown, ".png") |>
    png::readPNG() |>
    grid::rasterGrob(interpolate = TRUE) |>
    patchwork::wrap_elements()
  img_leftup <- paste0("../session1/stimuli/", img_leftup, ".png") |>
    png::readPNG() |>
    grid::rasterGrob(interpolate = TRUE) |>
    patchwork::wrap_elements()
  img_rightup <- paste0("../session1/stimuli/", img_rightup, ".png") |>
    png::readPNG() |>
    grid::rasterGrob(interpolate = TRUE) |>
    patchwork::wrap_elements()

  p <- ((p_err + theme(
    axis.title.x = element_blank(),
    plot.title = element_blank()
  )) /
    (p_rt + theme(plot.title = element_blank(), legend.position = "none"))) +
    patchwork::plot_layout(guides = "collect")

  wrap_elements(((img_leftup | patchwork::plot_spacer() | img_rightup) / p / (img_leftdown | patchwork::plot_spacer() | img_rightdown) +
    patchwork::plot_layout(heights = c(0.5, 1.5, 0.5)) +
    patchwork::plot_annotation(
      title = case_when(
        illname == "MullerLyer" ~ "Müller-Lyer",
        illname == "VerticalHorizontal" ~ "Vertical-Horizontal",
        TRUE ~ "Ebbinghaus"
      ),
      subtitle = question,
      theme = theme(plot.title = element_text(
        size=rel(1.75), face="bold", hjust=0.5),
        plot.subtitle = element_text(size=rel(1), face="italic", hjust=0.5, margin=margin(0,0,-30,0))))))
}

p_ebbinghaus <- plot_all(data, p_ebbinghaus_err, p_ebbinghaus_rt, "Which red circle is bigger?")
p_ebbinghaus

Müller-Lyer

Model Selection

test_models(filter(df, Illusion_Type == "MullerLyer"))

Error Rate

data <- filter(df, Illusion_Type == "MullerLyer")

Descriptive

plot_desc_errors(data)

Model Specification

formula <- brms::bf(
  Error ~ Illusion_Effect / (logmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))) +
    (1 + Illusion_Effect / (logmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))) | Participant),
  family = "bernoulli",
  decomp = "QR"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

illusion1_mullerlyer_err <- brms::brm(formula,
  data = data,
  refresh = 100,
  normalize = FALSE
)

save(illusion1_mullerlyer_err, file = "models/illusion1_mullerlyer_err.Rdata")

Model Inspection

load("models/illusion1_mullerlyer_err.Rdata")
load("models/gam_mullerlyer_err.Rdata")
p_mullerlyer_err <- plot_model_errors(data, illusion1_mullerlyer_err, gam_mullerlyer_err)
p_mullerlyer_err

Model Performance

performance::performance(illusion1_mullerlyer_err, metrics = c("R2", "ICC")) |>
  display()
R2 R2 (marg.) ICC
0.66 0.65 0.12

Response Time

data <- filter(df, Illusion_Type == "MullerLyer", Error == 0)

Descriptive

plot_desc_rt(data)

Model Specification

formula <- brms::bf(
  RT ~ Illusion_Effect / (sqrtmod(Illusion_Difference) * abs(Illusion_Strength)) + poly(ISI, 2) +
    (1 + Illusion_Effect / (sqrtmod(Illusion_Difference) * abs(Illusion_Strength)) | Participant),
  sigma ~ Illusion_Effect / (sqrtmod(Illusion_Difference) * abs(Illusion_Strength)) + poly(ISI, 2) +
    (1 | Participant),
  beta ~ Illusion_Effect / (sqrtmod(Illusion_Difference) * abs(Illusion_Strength)) + poly(ISI, 2) +
    (1 | Participant),
  family = "exgaussian",
  decomp = "QR"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

illusion1_mullerlyer_rt <- brms::brm(formula,
  data = data,
  refresh = 100,
  init = 0,
  normalize = FALSE
)

save(illusion1_mullerlyer_rt, file = "models/illusion1_mullerlyer_rt.Rdata")

Model Inspection

load("models/illusion1_mullerlyer_rt.Rdata")
load("models/gam_mullerlyer_rt.Rdata")
p_mullerlyer_rt <- plot_model_rt(data, illusion1_mullerlyer_rt, gam_mullerlyer_rt)
p_mullerlyer_rt

Model Performance

performance::performance(illusion1_mullerlyer_rt, metrics = c("R2", "ICC")) |>
  display()
R2 R2 (marg.)
0.47 0.14
plot_ppcheck(illusion1_mullerlyer_rt, gam_mullerlyer_rt)

Visualization

p_mullerlyer <- plot_all(data, p_mullerlyer_err, p_mullerlyer_rt, "Which red line is longer?")
p_mullerlyer

Vertical-Horizontal

Model Selection

test_models(filter(df, Illusion_Type == "VerticalHorizontal"))

Error Rate

data <- filter(df, Illusion_Type == "VerticalHorizontal")

Descriptive

plot_desc_errors(data)

Model Specification

formula <- brms::bf(
  Error ~ Illusion_Effect / (sqrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))) +
    (1 + Illusion_Effect / (sqrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))) | Participant),
  family = "bernoulli",
  decomp = "QR"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

illusion1_verticalhorizontal_err <- brms::brm(formula,
  data = data,
  refresh = 100,
  normalize = FALSE
)

save(illusion1_verticalhorizontal_err, file = "models/illusion1_verticalhorizontal_err.Rdata")

Model Inspection

load("models/illusion1_verticalhorizontal_err.Rdata")
load("models/gam_verticalhorizontal_err.Rdata")
p_verticalhorizontal_err <- plot_model_errors(
  data,
  illusion1_verticalhorizontal_err,
  gam_verticalhorizontal_err
)
p_verticalhorizontal_err

Model Performance

performance::performance(illusion1_verticalhorizontal_err, metrics = c("R2", "ICC")) |>
  display()
R2 R2 (marg.) ICC
0.44 0.37 0.06

Response Time

data <- filter(df, Illusion_Type == "VerticalHorizontal", Error == 0)

Descriptive

plot_desc_rt(data)

Model Specification

formula <- brms::bf(
  RT ~ Illusion_Effect / (cbrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))) + poly(ISI, 2) +
    (1 + Illusion_Effect / (cbrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))) | Participant),
  sigma ~ Illusion_Effect / (cbrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))) + poly(ISI, 2) +
    (1 | Participant),
  beta ~ Illusion_Effect / (cbrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))) + poly(ISI, 2) +
    (1 | Participant),
  family = "exgaussian",
  decomp = "QR"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

illusion1_verticalhorizontal_rt <- brms::brm(formula,
  data = data,
  refresh = 100,
  init = 0,
  normalize = FALSE
)

save(illusion1_verticalhorizontal_rt, file = "models/illusion1_verticalhorizontal_rt.Rdata")

Model Inspection

load("models/illusion1_verticalhorizontal_rt.Rdata")
load("models/gam_verticalhorizontal_rt.Rdata")
p_verticalhorizontal_rt <- plot_model_rt(
  data,
  illusion1_verticalhorizontal_rt,
  gam_verticalhorizontal_rt
)
p_verticalhorizontal_rt

Model Performance

performance::performance(illusion1_verticalhorizontal_rt, metrics = c("R2", "ICC")) |>
  display()
R2 R2 (marg.)
0.50 0.13
plot_ppcheck(illusion1_verticalhorizontal_rt, gam_verticalhorizontal_rt)

Visualization

p_verticalhorizontal <- plot_all(data, p_verticalhorizontal_err, p_verticalhorizontal_rt, "Which red line is longer?")
p_verticalhorizontal